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Abstract 

This study proposes a novel framework based on magnitude cumulant and surrogate analyses to 
reliably detect and estimate the intermittency coefficient from short-length coarse-resolution tur- 
bulent time series. Intermittency coefficients estimated from a large number of neutrally stratified 
atmospheric surface layer turbulent series from various field campaigns are shown to remarkably 
concur with well-known laboratory experimental results. In addition, surrogate-based hypothesis 
testing significantly reduces the likelihood of detecting a spurious non-zero intermittency coefficient 
from non-intermittent series. The discriminatory power of the proposed framework is promising for 
addressing the unresolved question of how atmospheric stability affects the intermittency properties 
of boundary layer turbulence. 
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I. INTRODUCTION 



Existence of small-scale intermittency is an intriguing yet unsettled topic in contemporary 
turbulence research. Over a number of decades, researchers had been trying to unravel inter- 
mittency in turbulence measurements and at the same time formulating diverse conceptual 
models to rationalize the observed intermittency.-^ Encouragingly, 'practical' implications 
of intermittency research outcomes are also being appreciated by the numerical turbulence 
modeling community and a critical knowledge transfer is taking place as evidenced by the 
recent literature.- - - 

One of the most widely used statistics characterizing the intermittent nature of turbulence 
is the so-called 'intermittency exponent' (/x).— From observational data, /i can be estimated 
directly or indirectly via several methods. The direct estimates typically involve appropriate 
characterization of the second-order scaling behavior of the local rate of energy dissipation 
(e) field. In this respect, several alternatives (e.g., second-order integral moment, two-point 
correlation function, spectral density) are available in the literature.- - — Recently, Cleve 
et al. 12 showed that among various direct approaches, the two-point correlation function 
(e(x + r)e(x)) of the energy dissipation field provides the most reliable estimates of \i. In 
this case, one can write 



where r is within the inertial range. Here, the angular brackets denote spatial averaging. 

The direct intermittency exponent estimation methods (based on e) require very high- 
resolution (resolving on the order of Kolmogorov scale) data series of pristine quality. Most 
commonly, fast-response hot wire measurements are used for this purpose.— _ — However, 
acquisition of hot wire data in a natural setting could be quite challenging. For example, 
in the case of atmospheric boundary layer (ABL) field experiments, one needs to perform 
meticulous hot wire calibration at short regular intervals in order to account for the ever 
changing (diurnally varying) ABL flow parameters.— 1 ^ 

The ABL community widely uses sonic anemometers for turbulent flux measurements 
in fields. In contrast to hot wires, these sensors require much lesser periodic calibration 
and maintenance. Unfortunately, path lengths (~ 10 cm) and sampling rate (~ 20 Hz) of 
conventional sonic anemometers are too coarse for direct /j estimation. In this paper, we 
will explore if sonic anemometer measurements, in lieu of hot wire data, can be reliably used 
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for indirect estimation of fi. 

One of the most popular indirect « estimation methods is associated with the scaling 
of sixth-order structure function.— liill With certain plausible assumptions, one can show 
that: 17 

3 6 (r) = ((u(x + r) - u{x)f) = ((Am) 6 ) ~ r 2 (e(x)e{x + r)) (2) 

Using Eqs. 1-2, one gets: ((Am) 6 ) ~ r 2 -M. Chambers and Antonia^ used this relatively 
simple indirect approach and obtained « ~ 0.2 in the atmospheric surface layer. From ade- 
quate statistical convergence standpoint, estimation of higher-order (specifically sixth-order 
in this context) structure functions require very long time series.— For instance, Cham- 
bers and AntoniaM used several runs of 15 min duration at a hot wire sampling frequency 
of ~ 1.2 kHz (i.e., ~ 1 million samples per series). On the other hand, 15-30 min sonic 
anemometer-based ABL turbulence series would typically consist of only 20,000 to 40,000 
samples. Understandably, the estimates of /i from sonic anemometer series using the tradi- 
tional sixth-order structure function approach will not be very reliable. This and some other 
theoretical considerations explained in the next section, motivate us to use an alternative es- 
timation approach, called the magnitude cumulant analysis, recently introduced by Delour et 
al.— . In this approach, only second-order magnitude cumulants (rather than data-intensive 
sixth-order structure functions) are needed to estimate the intermittency coefficient («). 
The inter-related objectives of this paper are twofold: 

(1) Assess the potential of the magnitude cumulant analysis in detecting and estimating 
intermittency from short-length coarse-resolution (sonic anemometer-acquired) ABL 
measurements; and 

(2) Design a rigorous hypothesis-testing framework which would reduce the likelihood of 
spurious detection of a non-zero « from non-intermittent (monofractal) series. 

The paper is structured as follows. In Section II, we briefly describe the magnitude 
cumulant analysis technique. The Iterative Amplitude Adjusted Fourier Transform (IAAFT) 
algorithm-based surrogate analysis, originally developed by the chaos theory community 
for detection of nonlinearity in time series,- 9 is shown to be very robust and reliable for 
intermittency hypothesis-testing. The IAAFT algorithm is presented in Section III. An 
extensive collection of observational data from various field campaigns is used in this study 
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and Section IV provides brief descriptions of these field datasets. Comprehensive results of 
intermittency estimation are presented in Section V and compared with published literature 
wherever possible. Lastly, Section VI summarizes our results and discusses perspectives for 
future research on the intriguing question of how atmospheric stability affects intermittency 
properties of boundary layer turbulence. 

II. MAGNITUDE CUMULANT ANALYSIS 

In the turbulence literature, the scaling exponent spectrum, ( q , is defined as: 



where S q (r) is the so-called g-th order structure function. As before, the angular bracket 
denotes spatial averaging and r is a separation distance that varies within the inertial range. 
If the scaling exponent ( q is a nonlinear function of q, then the field is called 'multifractal', 
otherwise it is termed 'monofractal'.- 1 ^ In the traditional structure function approach, esti- 
mation of « (= 2 — £ 6 ) requires a log- log plot of S 6 (r) vs r and subsequent extraction of the 
slope using a least-squares linear regression fit over a scaling regime (the inertial range). For 
short time series, computation of S& (r) is problematic due to statistical convergence. More- 
over, this problem is further compounded by the fact that even if the series is sufficiently long 
for statistical convergence of higher-order moments, the underlying nature of intermittency 
might theoretically limit the range over which the equivalency of statistical and geometrical 
interpretations of intermittency hold.— As a result, even accurate estimates of higher-order 
statistical moment will degenerate to a linear behavior of ( q for q larger than some q ma x 
prohibiting therefore an accurate estimation of intermittency using the structure function 
approach (see Lashermes et al.— for details). An alternative reliable method, first advocated 
by Delour et al.— , is to use the magnitude cumulant analysis. In this approach, the rela- 
tionship between the moments of velocity increments (Am) and the magnitude cumulants 
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(C„) reads as> 
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C 3 (r) = {(ln|Au|) 3 ) -3<(ln|A W |) 2 )<ln|Au|> 

(5c) 

+2(ln|Aw|) 3 ~ -c 3 ln(r) 
From Eqs. 3-5, it is straightforward to express the scaling exponent spectrum as: - 

°° q n 
n=l 

Furthermore, by invoking a relationship between velocity increments (Am) and local rate of 
dissipation fields (e) (similar to Eq. 2), one can arrive at>^ 

fi ~ 9c 2 (7) 

Therefore, estimation of the intermittency exponent \x would only require the computation 
of second-order magnitude cumulant, i.e., the second central moment of ln\Au\ (Eq. 5b). 

For large separation (r — > Li, where Li is the integral length scale), it is well documented 
that the probability distribution function (pdf) of velocity increments (Am) approaches a 
Gaussian distribution. For this scenario, the following results can be derived analytically: 

Ci(r) -> i(- 7 - ln(2)) = -0.64 (8a) 
C 2 (r) -> vr 2 /8 = 1.23 (8b) 
C 3 (r) - - 7 -m = -2.1 (8c) 

where 7 is the Euler Gamma constant = 0.577216, and £(3) is Apery's constant = 1.2020569. 
These asymptotic values of C\{r), Cz(r), and Cs(r) would be very useful to demarcate scaling 
regions in the case of short-length time series. 

It is noted that instead of a physical space-based magnitude cumulant analysis approach 
(i.e, Eqs. 4-5), one could also use wavelet-based magnitude cumulant analysis (see Venugopal 
et al.— for a geophysical application). A wavelet-based approach becomes necessary for 
nonstationary signals and signals with Holder exponents (h) outside the window of [0 1] 
In turbulence, (h) is close to K41 value of 1/3 and to best of our knowledge always found 
to be within the window of h £ [0 1].— & Thus, in the present study we decided to employ 
physical space-based magnitude cumulant analysis approach. 

Magnitude cumulant analysis of a synthetic fractional Brownian motion with h — 1/3 
(which displays K41 like k~ 5 ^ 3 spectrum) is shown in Fig. 1. The dashed line in C\(r) vs. 
ln(r) plot has the expected slope of 1/3. For almost the entire scaling range, both C 2 (r) and 



C 3 (r) remain close to the theoretical Gaussian values of 1.23 and —2.1, respectively. This 
signal does not show any sign of multifractality (expected) as the slope of C 2 (r) vs. ln(r) 
cannot be claimed to be different from zero. 




-1 0123456 -1 0123456 

ln(r) ln(r) 

FIG. 1: A synthetic fractional Brownian motion (h = 1/3) (top-left). The magnitude cumulants 
C\{r) (top-right), G 2 (r) (bottom- left), and Cs(r) (bottom-right) are also shown. The dashed line 
in the top-right subplot shows the slope —c\ = 1/3. 



III. SURROGATE ANALYSIS 

Noise is omnipresent in any measured signal and turbulence signals are no exceptions. 
In addition to noise, limited amount of data (finite sample settings) in most field mea- 
surements could challenge intermittency detection and estimation even with the magnitude 
cumulant analysis method (e.g., assessment of a small non-zero slope in the G 2 (r) vs. ln(r) 
plots). In this paper, we utilize a hypothesis-testing framework, based on surrogate analysis, 
in conjunction with magnitude cumulant analysis, for detecting and accurately estimating 
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intermittency from short-length sonic anemometer measurements. 

The concept of surrogates (stochastic realizations which preserve only certain character- 
istics of a process) was introduced into the chaos theory literature to provide a rigorous 
statistical test for the null hypothesis that an observed time series has been generated by 
a linear stochastic process (see Theiler et al.,— Kantz and Schreiber,— Basu and Foufoula- 
Georgiou,— and the references there in). Over the years, several varieties of surrogates 
(randomly shuffled surrogates, Fourier phase randomized surrogates, iterative amplitude 
adjusted Fourier transform - IAAFT surrogates, stochastic IAAFT surrogates, and so on) 
have been proposed in the literature.— In this paper, we will use the IAAFT algorithm pro- 
posed by Schreiber and Schmitz. 19 IAAFT surrogates preserve the correlation structure (thus 
power spectrum owing to Wiener-Khinchin theorem) and the probability density function 
of a given time series. Apart from nonlinearity detection,—^ the IAAFT surrogates have 
also been used to define a precipitation forecast quality index,— and to generate synthetic 
cloud fields.— 

In the turbulence literature, surrogate analysis-based hypothesis-testing is virtually 
nonexistent with an exception of the paper by Nikora et al.— They used simple Fourier 
phase randomization approach (pdf of the original turbulence series was not preserved) in 
identifying the effects of turbulence intermittency and spectral energy flux. In comparison to 
the Fourier phase randomized approach, the IAAFT algorithm used in the present study de- 
signs stronger statistical test (owing to its ability to preserve the integral pdf of the original 
signal) for the null hypothesis that an observed turbulence series is non- intermittent. 

In Fig. 2, a sonic anemometer turbulence series and its IAAFT surrogate are shown. 
By construction, they have the same pdf (bottom-right plot of Fig. 2) and virtually indis- 
tinguishable autocorrelation function (bottom-left plot of Fig. 2). Basic properties of the 
original turbulence series and its surrogate are provided in Table I. Tj and Lj denote integral 
time and length scales, respectively. 



where R(t) is the autocorrelation function. From Fig. 2 and Table I, we can safely infer that 
the IAAFT surrogate captures the integral pdf and autocorrelation function of the original 
velocity series rather accurately. Later on, we will show that the IAAFT surrogates do not 




(9a) 
(9b) 
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have the ability to capture the scale-dependent pdfs of velocity increments and this forms 
the basis for the proposed intermittency hypothesis-testing. 




FIG. 2: Longitudinal velocity time series measured by a sonic anemometer (top-left) and its surro- 
gate series generated by IAAFT methodology (top-right). Both series have approximately the same 
autocorrelation (bottom- left) and exactly the same probability density function (bottom-right). 



TABLE I: Basic Statistics of a Sonic Anemometer Turbulence Series and its Surrogate 



Series Type U (ms 1 ) 


a u (ms" 1 ) 


Ti (s) 


Li (m) 


Turbulence 3.73 


1.06 


1.47 


5.47 


Surrogate 3.73 


1.06 


1.46 


5.43 
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IV. DESCRIPTION OF DATA 



In this study, we primarily made use of an extensive atmospheric boundary layer tur- 
bulence dataset (comprising of sonic anemometer measurements) collected by various re- 
searchers from Johns Hopkins University, and the University of California-Davis during 
Davis 1994, 1995, 1996, 1999 field studies. Comprehensive description of these field ex- 
periments (e.g., surface cover, fetch, instrumentation, sampling frequency) can be found in 
Pahlow et al. 37 Briefly, the collective attributes of the field dataset explored in this study 
are as follows: (i) surface cover: bare soil, and beans; (ii) sampling frequency: 18 to 21 Hz; 
(iii) sampling period: 20 to 30 minutes; (iv) sensor height (z): 0.96 to 4.28 m. 

The ABL field measurements are seldom free from mesoscale disturbances, wave activities, 
nonstationarities etc. The situation could be further aggravated by several kinds of sensor 
errors (e.g., random spikes, amplitude resolution error, drop outs, discontinuities etc.). Thus, 
stringent quality control and preprocessing of field data is of utmost importance for any 
rigorous statistical analysis. Our quality control and preprocessing strategies are described 
in detail in Basu et al.— After the quality control and preprocessing steps, we were left with 
139 'reliable' near-neutral (\z/L\ < 0.05, where z is the sensor height and L denotes the 
Monin-Obukhov length) sets of runs for estimating the intermittency exponents. 

We also estimated fi from a fast-response (10 kHz) hot wire ABL turbulence series utilizing 
the magnitude cumulant analysis. The hot wire measurements were taken at the Surface 
Layer Turbulence and Environmental Science Test (SLTEST) facility located in the western 
Utah Great Salt Lake desert under near-neutral atmospheric condition.—^ In the following 
section, we will show that the intermittency exponent and other relevant statistics derived 
from this high Reynolds number (-Re) hot wire measurement are surprisingly similar to 
various published lower Re laboratory experimental findings, and serve as benchmarks in 
the present study. 



TABLE II: Mean Flow Characteristics of the Field Measurements 



Sensor Type 


z (m) 


U (ms 1 ) a u (ms x ) 


Ti (s) 


U 


(m) 


Hot wire Anemometer 


2.01 


5.99 0.74 


5.71 


34.22 


Sonic Anemometer 


0.96 - 4.28 1.60 - 7.30 0.34 - 1.57 1.08 - 9.09 


2.57- 


- 33.22 
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Mean flow characteristics of all the field measurements are given in Table II. For all the 
analyses, we have invoked Taylor's hypothesis to convert time series to spatial series. 

V. RESULTS 

A. Analysis of Hot Wire Measurements 

In this section, hot wire measurements and their surrogates are analyzed to: (a) demon- 
strate the ability of the magnitude cumulant analysis to accurately estimate the intermit- 
tency structure of turbulent velocity series, and (b) establish that the surrogate series, while 
preserving the pdf and spectrum of the original data, destroy the intermittency structure. 

In Fig. 3, the magnitude cumulants computed from the hot wire measurements of Kunkel 
and Marusic - are shown. This turbulence series is 30 min long (~ 18 million data points) 
and captures scales down to the Kolmogorov scale. The following observations can be made 
from Fig. 3: 

• c t ) urb computed from the turbulence velocity signal is close to —0.36. It agrees quite 
well with the existing results from low Re laboratory experiments.— 1 ^ 

• cty rb is approximately equal to 0.03. Delour et al.— and Chevillard et al.— reported 
c 2 = 0.025 ± 0.003 based on several experiments and claimed it to be 'universal'. 
From Eq. 7, we can compute the intermittency exponent /i ~ 9 • c| urb ~ 0.27. In 
the literature, researchers have reported /x ranging from 0.18 to 0.7.—^ In the case of 
atmospheric data, the 'best' direct estimate is 0.25±0.05^ and our indirect magnitude 
cumulant analysis-based result is in agreement with it. 

• From Eq. 6, we can derive: C 2 urb = -2 ■ c\ urb - 2 ■ c l 2 urb = 2 ■ 0.36 - 2 ■ 0.03 = 0.66. Fig. 3 
(bottom-right) shows the second-order structure function. The slope of this plot gives 
^turb _ q qq Thus, our results are self-consistent. 

• In the inertial range, c t ^ ,rb seems to be zero. This indicates that the statistics of the 
velocity increments are possibly log-normal.— 1 ^ 
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FIG. 3: C\(r/Li) (top-left), C2(r/Lj) (top-right) and C^{r/Li) (bottom-left) computed using the 
hot wire measurements of Kunkel and Marusic.- 6 The second-order structure function is also shown 
(bottom-right). Extended range of inertial-range is evident in all the sub-plots. The circles denote 
the statistics corresponding to the original turbulence series and the stars represent the statistics 
computed from the IAAFT surrogate series. Clearly, the original series portray the signatures of 
a multifractal process. In contrast, the surrogate series shows the signs of monofractality. The 
dashed lines show the slopes -c[ urb = 0.36 (top-left), -cf rh = -0.03 (top-right), and C^ rh = 0.66 
(bottom-right), respectively. 

• c S 2 TT estimated from the surrogate of the measured turbulence velocity series is zero, 
i.e., the surrogate series is non-intermittent. 

• By construction, the surrogate series (i.e., turbulence without intermittency) preserves 
the second-order statistics. So, Cj^ rb = Q urr . Using this relationship, the fact that 
c surr _ anc j j£)q_ g ; j s straightforward to show that: c t ^ rh + c 2 urfe = c\ urr . In the 
present case this rsults in: cf rr = —0.36 + 0.03 = —0.33 (see also Fig. 3). Thus, 
Cr rr = -c\ urr = 0.33 is in full accord with K41 hypothesis of Ci = 1/3. 1 
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The relationship c l ^ rb + c% irb = c{ urr has significant practical implication. It insinuates 
that one can roughly estimate c^" -6 (and thus //) by means of first-order magnitude 
cumulants of turbulence and its corresponding surrogate, i.e., d^ rh = d\ rb — c{ urr . 
In our opinion, for estimating intermittency in short length geophysical signals, this 
simple indirect method which does not require even second-order magnitude cumulant 
computation would be quite useful. 

• In the turbulence literature, there is a general consensus that £3 = 1. From our 
results, we find C^ urb = -3 • c[ urb - 9/2 • c^ rb = 1.08 - 0.135 = 0.95, very close to the 
well-accepted value. 

We proceed further by comparing the pdf of the velocity and surrogate increments using 
the skewness, asymmetry factor, and flatness defined as: 

((Am) 3 ) 

Skewness(r) = ... , 9 ,„ /9 (10a) 

((Am) 3 ) , , 

Asymmetry [r) = TTT^T^K (10b) 

Flatness(r) = (10c) 

From Fig. 4 (left), it is evident that the original turbulence increment series show negative 
skewness (up to ~ 0.6) for small scales in accord with the existing literature (e.g., Chevillard 
et al.—). This negative skewness is believed to be related to the vortex folding and stretching 
process. Malecot et al.— argued that the asymmetry factor (see Eq. 10 for definition) is 
a better measure of the asymmetry of the pdf than the skewness. We found that both 
of these signed odd-order moments behave quite similarly (Fig. 4, left). The origin of 
spurious oscillations of these odd-order moments for large scales (ln{r/Li) > —4) is not well 
understood. The flatness plot (Figure 4, right) also portrays anticipated characteristics. 
Flatness corresponding to the integral scale is close to 3 (hallmark of Gaussian velocity 
increments) and becomes exceedingly large for smaller scales. In contrast, the surrogate 
shows Gaussian characteristics for all scales. This corroborates the fact that surrogates 
cannot capture the pdfs of turbulence velocity increments. 
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FIG. 4: Negative skewness, asymmetry factor (left), and flatness (right) of the longitudinal velocity 
increments and the increments of the surrogate series. We utilized the hot wire measurements of 
Kunkel and Marusic. 16 

B. Analysis of Sonic Anemometer Measurements 

Magnitude cumulants and second-order structure functions computed from a sonic 
anemometer series are shown in Fig. 5. The trends are very similar to Fig. 3, albeit 
quite noisy. From this figure, we calculated C\ = —0.37, c 2 = 0.06, and (2 = 0.63. We 
would like to emphasize that even in this short time series scenario, we can reliably detect 
intermittency with the help of IAAFT surrogate (see Fig. 5 top-right). Admittedly, the 
estimation of C2 is possibly not very accurate. The estimation can be improved by using 
quenched or annealed averaging strategy (discussed below). 



It is quite difficult to manually yet objectively select scaling ranges from a large dataset 
(specifically 139 near- neutral turbulence series). So, we used an automated scaling range 
of [4 • U/f s Li/2] for individual series. Here U, f s , and Lj denote mean velocity, sampling 
frequency, and integral length scale, respectively. From each individual turbulence series and 
their surrogates, we calculated the corresponding c\ and c 2 values. Subsequently, from these 
139 c* , c| ', c{ urr , and c s 2 urr combinations, we computed the best estimates (ensemble 
mean plus/minus one standard deviation) as: (c' urfe ) = — 0.33 =F 0.03, (c t 2 urb ) = 0.038 ±0.017, 
(cl urr ) = -0.30 =p 0.03, and (c s 2 urr ) = 0.002 ± 0.013. This averaging strategy is similar to 
the quenched averaging method used by Arneodo et al.— The key result: (c t 2 urb ) >> (c 2 urr ), 
without any doubt, once again guarantees that the IAAFT surrogates can be faithfully used 
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FIG. 5: Same as Figure 3, but estimated from a sonic anemometer series. 

for turbulence intermittency detection testing. 

As an alternative strategy, using the annealed averaging method-^, we have also com- 
puted the average of the magnitude cumulants (i.e., (Ci (r/Lj)), (C 2 (r/Li)) ) from the 
same turbulence and surrogate datasets (Fig. 6). From this figure, we estimate the slopes 
as: c\ urb = -0.35, c l 2 urb = 0.042, c\ urr = -0.31, c s 2 urr = 0.002. Obviously, there is no sig- 
nificant discrepancy between the quenched and annealed averaged statistics, as would be 
wishfully expected. Lastly, these statistics highlight that the relationship established in 
Section V-A, i.e., c' Mrfc + c 2 urb = c{ urr , is also valid under annealed averaging. 



VI. CONCLUDING REMARKS AND FUTURE PERSPECTIVES 

In this work, we have established a framework based on magnitude cumulant and surro- 
gate analyses to reliably detect and estimate the intermittency coefficient from short turbu- 
lent time series. By virtue of this framework, ensemble scaling results extracted from a large 
number of neutrally stratified atmospheric surface layer turbulent series (predominantly ac- 
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FIG. 6: Ensemble averaged Ci(r/Lj) and (^(r/Lj) plots from 139 near-neutral turbulence series 
(circles) and corresponding surrogates (stars). The dashed lines show the slopes — c^ rfe = 0.35 
(left), and -c^ = -0.042 (right), respectively. 

quired by slow response sonic anemometers) from various field campaigns remarkably concur 
with well-known published (mostly laboratory experimental) results. 

The focus of the present study was on neutrally stratified atmospheric turbulence. How- 
ever, it is widely known that neutral stability conditions are rarely encountered in the real 
atmosphere. Most of the time, the atmospheric boundary layer is strongly modulated by 
buoyancy. It is commonly assumed that the effect of atmospheric stability is felt only at 
the 'buoyancy range', which has scales considerably larger than the inertial range. Recently, 
Aivalis et al.— studied the intermittency behavior of temperature in the convective surface 
layer using cold wire anemometry. They found that the classical inertial range remains in- 
tact in convective surface layer and the scaling exponents approach values appropriate to 
the intermittent case of isotropic turbulence. They also noticed that the scaling exponents 
corresponding to the buoyancy range are highly anomalous. In contrast, Shi et al.— found 
that in the case of temperature, the inertial range scaling exponents are unambiguously 
impacted by atmospheric stability. However, their results are inconclusive in the case of 
velocity signals under different stability regimes. We believe that the unresolved question 
as to whether or not the inertial-range intermittency is influenced by large-scale anisotropic 
forcing of atmospheric stability is of great consequence and needs further consideration. 
Abundant high-quality slow response data collected under different atmospheric regimes in 
recent years (e.g., Cooperative Atmosphere-Surface Exchange Study - CASES99 42 ) could be 
coupled with the robust scaling analysis and estimation framework explored in this study in 
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order to shed new light into this fundamental problem. 
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